BuildNetwork

library(AnNet)
library(synaptome.db)
library(pander)
library(ggplot2)
library(ggrepel)

scale <- function(x, VALUE=NULL){

    x = as.numeric(as.vector(x))

    xmin <- min(x,na.rm=T)
    xmax <- max(x,na.rm=T)
   
    if( is.null(VALUE) ){

        x  <- x-xmin
        x  <- ifelse(!is.na(x), x/(xmax-xmin), NA) 

        return(x)
    }

    value = as.numeric(as.vector(value)[1])
    value = value-xmin
    value = ifelse(!is.na(value), value/(xmax-xmin), NA) 
    return(value)
}

Title AnNet: library for comprehensive network analysis

Authors: Colin Mclean, Anatoly Sorokin, J. Douglas Armstrong and Oksana Sorokina

#Introduction

Overview of capabilities

Build the network

The package allows building the network from the data frame, where the rows correspond to the edges of the graph, or for the list of nodes (genes), for which the information will be retrieved from the SynaptomeDB.

Build network for the node list extracted from SynaptomeDB


cid<-match('Presynaptic',getCompartments()$Name) # Let's get the ID for presynaptic compartment
cid
#> [1] 2
t<-getAllGenes4Compartment(cid) # Now we need to collect all the gene IDs for presinaptic  compartment
dim(t)
#> [1] 2304    8
head(t)
#> # A tibble: 6 × 8
#>   GeneID MGI       HumanEntrez MouseEntrez RatEntrez HumanName MouseName RatName
#>    <int> <chr>           <int>       <int>     <int> <chr>     <chr>     <chr>  
#> 1      1 MGI:1277…        1742       13385     29495 DLG4      Dlg4      Dlg4   
#> 2      2 MGI:88256         815       12322     25400 CAMK2A    Camk2a    Camk2a 
#> 3      3 MGI:96568        9118      226180     24503 INA       Ina       Ina    
#> 4      4 MGI:98388        6711       20742    305614 SPTBN1    Sptbn1    Sptbn1 
#> 5      5 MGI:88257         816       12323     24245 CAMK2B    Camk2b    Camk2b 
#> 6      6 MGI:1344…        1740       23859     64053 DLG2      Dlg2      Dlg2
gg<-buildFromSynaptomeByEntrez(t$HumanEntrez) # finally, build the graph using respecctive gene EntrezIDs as node IDs
summary(gg)
#> IGRAPH be669ac UN-- 1780 6620 -- 
#> + attr: name (v/c)

Annotate the nodes with node attributes

Gene name


gg<-annotateGeneNames(gg)
#> 'select()' returned 1:1 mapping between keys and columns
summary(gg)
#> IGRAPH be669ac UN-- 1780 6620 -- 
#> + attr: name (v/c), GeneName (v/c)
head(V(gg))
#> + 6/1780 vertices, named, from be669ac:
#> [1] 273   6455  1759  1785  26052 2923
head(V(gg)$GeneName)
#> [1] "AMPH"   "SH3GL1" "DNM1"   "DNM2"   "DNM3"   "PDIA3"

Clustering

alg<- 'louvain'
gg<-calcClustering(gg,alg = alg)
{r cluster.graph, cache=TRUE}
gg<-calcAllClustering(gg)

Centrality

gg <- calcCentrality(gg)

Consensus matrix

Build consensus matrix for louvain clustering

conmat<-makeConsensusMatrix(gg,N=100,mask = 10,alg = alg,type = 2)
steps <- 100
Fn  <- ecdf(conmat[lower.tri(conmat)])
X<-seq(0,1,length.out=steps+1)
cdf<-Fn(X)
dt<-data.frame(cons=X,cdf=cdf)
ggplot(dt,aes(x=cons,y=cdf))+geom_line()+
      theme(            
        axis.title.x=element_text(face="bold",size=rel(2.5)),
        axis.title.y=element_text(face="bold",size=rel(2.5)),
        legend.title=element_text(face="bold",size=rel(1.5)),
        legend.text=element_text(face="bold",size=rel(1.5)),
        legend.key=element_blank())+
    theme(panel.grid.major = element_line(colour="grey40",size=0.2),
          panel.grid.minor = element_line(colour="grey40",size=0.1),
          panel.background = element_rect(fill="white"),
          panel.border = element_rect(linetype="solid",fill=NA))

Cluster robustness

clrob<-getRobustness(gg,alg,conmat)
pander(clrob)
alg C Cn Crob CrobScaled
louvain 1 189 0.09604 1
louvain 2 217 0.09 0.4831
louvain 3 282 0.08741 0.2615
louvain 4 139 0.08849 0.3533
louvain 5 121 0.09316 0.7536
louvain 6 163 0.08759 0.2769
louvain 7 30 0.08436 0
louvain 8 225 0.08831 0.338
louvain 9 150 0.0866 0.1922
louvain 10 108 0.09031 0.5097
louvain 11 106 0.08519 0.0714
louvain 12 19 0.09012 0.4934
louvain 13 31 0.09251 0.6978

Bridgeness

bm<-getBridgeness(gg,alg,conmat)
#> calculating Bridgeness for:  louvain
pander(head(bm))
ENTREZ.ID GENE.NAME BRIDGENESS.louvain
273 AMPH 0.0957579145556555
6455 SH3GL1 0.452106501493934
1759 DNM1 0.0892319279050531
1785 DNM2 0.254828113224788
26052 DNM3 0.147644688926001
2923 PDIA3 0.548658380347371 
VIPs=c('8495','22999','8927','8573','26059','8497','27445','8499')
# VIPs=c('81876','10890','51552','5874','5862','11021','54734','5865','5864',
#            '9522','192683','10067','10396','9296','527','9114','537','535',
#            '528','51382','534','51606','523','80331','114569','127262','57084',
#            '57030','388662','6853','6854','8224','9900','9899','9145','9143',
#            '6855','132204','6857','127833','6861','529','526','140679','7781',
#            '81615','6844','6843')
indx   <- match(V(gg)$name,VIPs)
group <- ifelse( is.na(indx), 0,1)
MainDivSize <- 0.8
xmin        <- 0
xmax        <- 1
ymin        <- 0
ymax        <- 1
Xlab <- "Semilocal Centrality (SL)" 
Ylab <- "Bridgeness (B)"
X    <- as.numeric(igraph::get.vertex.attribute(gg,"SL",V(gg)))
X    <- scale(X)
Y       <- as.numeric(as.vector(bm[,dim(bm)[2]])) 
lbls <- ifelse(!is.na(indx),V(gg)$GeneName,"")
dt<-data.frame(X=X,Y=Y,vips=group,entres=V(gg)$name,name=V(gg)$GeneName)
dt_vips<-dt[dt$vips==1,]
dt_res<-dt[dt$vips==0,]
##--- baseColor of genes
baseColor="royalblue2"

##--- SPcolor, colour highlighting any 'specical' genes
SPColor="royalblue2"

##--- PSDColor, colour of core PSD95 interactor genes
PSDColor="magenta"

ggplot(dt,aes(x=X,y=Y,label=name))+#geom_point()+
    geom_point(data=dt_vips,
               aes(x=X,y=Y),colour=baseColor,size = 7,shape=15,show.legend=F)+
    geom_point(data=dt_res,
               aes(x=X,y=Y, alpha=(X*Y)), size = 3,shape=16,show.legend=F)+
    geom_label_repel(aes(label=as.vector(lbls)),
                     fontface='bold',color='black',fill='white',box.padding=0.1,
                     point.padding=NA,label.padding=0.15,segment.color='black',
                     force=1,size=rel(3.8),show.legend=F,max.overlaps=200)+
      labs(x=Xlab,y=Ylab,title=sprintf("%s",alg))+
    scale_x_continuous(expand = c(0, 0), limits = c(xmin, xmax)) + 
    scale_y_continuous(expand = c(0, 0), limits = c(ymin, ymax))+
    theme(            
        axis.title.x=element_text(face="bold",size=rel(2.5)),
        axis.title.y=element_text(face="bold",size=rel(2.5)),
        legend.title=element_text(face="bold",size=rel(1.5)),
        legend.text=element_text(face="bold",size=rel(1.5)),
        legend.key=element_blank())+
    theme(panel.grid.major = element_line(colour="grey40",size=0.2),
          panel.grid.minor = element_line(colour="grey40",size=0.1),
          panel.background = element_rect(fill="white"),
          panel.border = element_rect(linetype="solid",fill=NA))+
        geom_vline(xintercept=0.5,colour="grey40",size=MainDivSize,linetype=2,show.legend=F)+
    geom_hline(yintercept=0.5,colour="grey40",size=MainDivSize,linetype=2,show.legend=F)
#> Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
#> increasing max.overlaps

library(plotly)
p<-ggplot(dt,aes(x=X,y=Y,label=name))+#geom_point()+
    geom_point(data=dt_vips,
               aes(x=X,y=Y),colour=baseColor,shape=15,show.legend=F)+
    geom_point(data=dt_res,
               aes(x=X,y=Y, alpha=(X*Y)),shape=16,show.legend=F)+
      labs(x=Xlab,y=Ylab,title=sprintf("%s",alg))+
    scale_x_continuous(expand = c(0, 0), limits = c(xmin, xmax)) + 
    scale_y_continuous(expand = c(0, 0), limits = c(ymin, ymax))+
    theme(            
        axis.title.x=element_text(face="bold",size=rel(2.5)),
        axis.title.y=element_text(face="bold",size=rel(2.5)),
        legend.title=element_text(face="bold",size=rel(1.5)),
        legend.text=element_text(face="bold",size=rel(1.5)),
        legend.key=element_blank())+
    theme(panel.grid.major = element_line(colour="grey40",size=0.2),
          panel.grid.minor = element_line(colour="grey40",size=0.1),
          panel.background = element_rect(fill="white"),
          panel.border = element_rect(linetype="solid",fill=NA))+
        geom_vline(xintercept=0.5,colour="grey40",size=MainDivSize,linetype=2,show.legend=F)+
    geom_hline(yintercept=0.5,colour="grey40",size=MainDivSize,linetype=2,show.legend=F)
ggplotly(p)
#> Warning: `gather_()` was deprecated in tidyr 1.2.0.
#> Please use `gather()` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.